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Abstract . The well-known Waddington's epigenetic landscape of cell-fate determination is 
not static but varies because of the dynamic gene regulation during development. However, 
existing mathematical models with few state variables and fixed parameters are inadequate in 
characterizing the temporal transformation of the landscape. Here we simulate a decision- 
switch model of gene regulation with more than two state variables and with time-varying 
repression among regulatory factors. We are able to demonstrate multi-lineage differentiation 
at different timescales that portrays the branching canals in Waddington's illustration. We 
also present a repressilator-type system that activates suppressed genes via sustained 
oscillations in a flattened landscape, hence providing an alternative scheme for cellular 
reprogramming. The time-dependent parameters governed by gradient-based dynamics 
regulate cell differentiation, dedifferentiation and transdifferentiation. Our prediction 
assimilates the theories of branching and structural oscillations in cell-fate determination, 
which reveals a potential blueprint for cell differentiation and associated diseases, such as 
cancer. 

Keywords , gene regulatory network, stem cells, pluripotency, synthetic biology, 
multistability, attractor 
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Main text: 

Waddington's epigenetic landscape illustrates the canalization in the cell differentiation and 
fate determination process 1 " 4 . The topography of Waddington's illustration represents the 
developmental pathways of tissues formed from totipotent and pluripotent cells to terminally- 
differentiated specialized cells (Fig. la). Various theoretical studies have quantified 
Waddington's epigenetic landscape and are able to predict bistability in gene regulatory 
networks (GRNs) 5 " 9 . However, many of the mathematical models only consider at most two 
regulatory factors, and focus on static epigenetic landscape represented by fixed parameter 
values. In reality, the topography of Waddington's illustration is dynamic and high- 
dimensional, and the parameters that represent gene regulation are changing during the 
development of an organism 10 " 16 . Mathematical models with two regulatory factors and fixed 
parameter values only describe a particular temporal scenario in cell differentiation. 

The mechanisms that regulate gene expression, such as kinetics of gene regulatory 
factors (GRFs) and the structure of GRF-GRF interaction, influence the outcome of cell-fate 
determination 7 ' 10 ' 11 ' 14 ' 15 ' 17 . Waddington observed that changes in these mechanisms could alter 
the epigenetic landscape leading to cell-lineage switching 1 . The changes in the GRF-GRF 
interaction do not necessarily entail mutations but can be due to normal processes. 
Mathematically, the variations in gene regulation can be represented by modifications in the 
parameter values of the quantitative models. Bifurcation analyses of existing models have 
been done 5 ' 7 ' 11 ' 15 ' 16 , but most of them do not provide elaborate illustrations of cells trailing the 
high-dimensional dynamic pathways. Here we present numerical illustrations of cells trailing 
different epigenetic routes such that the pathways transform due to changes in the strength of 
repressive interaction among multiple GRFs (see Box 1 and Fig. lb for the mathematical 
model). The GRFs in the model (Fig. lb) has mutual repression because a mature cell 
expresses only one phenotype and constrains the expression of the other phenotypes. 
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A desired cell fate can be a cell type/phenotype that is essential for proper normal 
development, or desired cell type during cellular engineering. Our main assumption is that 
gradient-based optimization governs the transformations of the pathways leading to the 
desired cell fate, following the theory that differentiating cells choose the steeper pathways 
(canals) in the epigenetic landscape. This assumption assures that the cells trail the nearby 
steepest pathways shaped by the time-varying antagonistic interaction among the GRFs (see 
Methods). The gradient-based method can be considered as a cell-fate induction strategy such 
that the cells move from a pluripotent state, which has higher entropy, towards differentiated 

1 8 

state with lesser entropy . We show that the dynamic GRF-GRF interaction can illustrate the 
cascade of branching canals in Waddington's illustration. It can also describe cell plasticity 
by allowing cell-lineage switching, such as by transdifferentiation and dedifferentiation. 



Box 1: Mathematical model 

We consider a minimal gene regulatory network (GRN) of the form shown in Figure lb, 
where the number of nodes (n) is arbitrary. This GRN is a representation of the network that 

10 22 27 34 

characterizes decision switches in cell-fate determination ' ' ' . For simplicity, we draw the 
network in such a way that each node represents a gene regulatory factor (GRF) involved in 
expressing a specific cell type/phenotype 10 ' 12 . A node in this GRN represents either a specific 
GRF, or a coarse-grained subnetwork of multiple factors that can be treated as single GRF. 
An example of a GRN in the form shown in Fig. lb is the coarse-grained mesenchymal 
transcription network with RUNX2, PPAR-y and SOX9 as nodes 17 . 

One of the simple high-dimensional models that describe the qualitative dynamics of 
the GRN is the following system of differential equations: 

^ = ^ j-pX i + g,i = l,2,...,n. (BoxEq. 1) 
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The state variable X, represents the strength/concentration of the GRF involved in expressing 
the gene towards the i-th cell type. The parameters ft>0, p>0, g>0 and y,y>0 are the efficiency 
of GRF in expressing the corresponding gene, degradation rate, basal constitutive growth rate, 
and time-varying interaction coefficient associated with the inhibition of X, by Xj, 
respectively. The kinetics of GRF auto-activation is a sigmoidal increasing function which is 
negatively influenced by the strength of the other GRFs. This model is originally proposed by 
Cinquin and Demongeot 10 ' 12 but their model has symmetric parameters unlike Box Eq. 1 
which admits asymmetric repression coefficient 

We assume that the goal of gene regulation is to maximize the strength of GRFs so 
that the outcome is moving towards the steepest canal and possibly towards the deepest 
valley in the epigenetic landscape, which is expected in cell-fate determination. Here we 
apply a gradient-based optimization method to the time-varying repression strength y 2 y to 
drive cell-fate induction in the direction of the nearby steep canal (see Methods). This method 
is governed by a timescale factor that declines through time 5 ' 11 . 

The equilibrium points and sustained oscillations of this model lie on the hyperspace 



8 8 + P 



P P J 
(End of Box 1) 



The model has at most 3 n equilibrium points 



20 



Results . We observe two significant dynamics in our simulations. The first one is multi- 
lineage differentiation via sequentially branching developmental pathways. This sequential 
branching portrays the canalization in Waddington's landscape at different timescales. The 
pathways trailed by the differentiating cells depend not only on the structure of the GRN and 
parameter values but also on the initial condition (see Supplementary Fig. 1). The second one 
is flattening of the epigenetic landscape which eliminates the deep valleys, resulting in 
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sustained oscillations. The GRN that generates the oscillations is an attracting oscillator 
where it can lift the strength of a suppressed GRF. This oscillatory behavior is structural in 
nature (i.e., based on the structure of the GRN) and not due to noise nor delay. Therefore, we 
are able to integrate into a single theory the theories of branching and structural oscillations 
in cell-fate determination. 

Mathematically, the sequential branching in the epigenetic landscape towards 
different cell types represents convergence to one of the equilibrium points (Supplementary 
Fig. 1). The variations in the interaction coefficients drive changes in the topography of the 
landscape (Supplementary Fig. 2) and eventually stabilize at different timescales, hence 
sequential branching become possible (refer to the timescale factor in Methods). The pathway 
bifurcates from one branch (a primed state) to multiple branches, which may further have 
sub-branches that culminate in branch endpoints (Figs. 2a, 2c-2e and Supplementary Fig. 21). 
The final structure of the GRN, which is determined by the parameters stabilized at different 
timescales, dictates the number and location of the branch endpoints (Figs. 2a-2e; see 
Supplementary Fig. 21 for a system with 10 GRFs). Generally, in order for sequential 
branching to arise, the initial structure of the GRN and parameter values should allow 
bistability or multistability which is a property of multipotent and pluripotent cells 
(Supplementary Fig. 1). The entropy and the quasi-potential of the landscape are lower in 
differentiated cells compared to the undifferentiated state as expected (e.g., Fig. 2f and see 
the quasi-potential axis in Figs. 2a and 2b). In addition, the pathways with two endpoints that 
are distant from each other are usually more robust against stochastic noise compared to the 
pathways with many endpoints (Supplementary Figs. 3, 6, 8 and 10). This implies that cells 
trailing the pathways with more endpoints can be candidates for stochastic direct 
transdifferentiation . 
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There are cases where a stable equilibrium point vanishes and a stable limit cycle 
emerges, especially when the repressive interaction among GRFs is asymmetric. This stable 
limit cycle is generated by an oscillator that attracts suppressed genes in partially or 
terminally-differentiated cells, resulting in activation with fluctuating kinetics (Figs. 3b and 
3c). One example of an oscillator is a repressilator-type network (Fig. 3a), which is similar 
but not exactly identical to the repressilator proposed by Elowitz and Leibler 19 . In this 
repressilator-type network, the strength of repression in one loop is stronger than the strength 
of repression in the reverse loop (asymmetric interaction). The sustained oscillations 
generated by this repressilator-type network can arise in a GRN with three or more nodes 
(odd or even number of GRFs; e.g., Figs. 3b and 3c). The dynamics of this oscillator can be 
illustrated by an epigenetic landscape with flattened topography, that is, there are no deep 
valleys in the route of the differentiating cells, and the cells are continually sliding in zigzag 
canals without endpoint. 

The oscillations drive the strengths of the GRFs to have alternating positive and 
negative rates, whereas the gradient-based optimization method forces the dynamics towards 
positive rates only. Thus, oscillatory behavior is not optimal in the sense of differentiation 
towards cell types located at deep valleys in the landscape. We then expect that gradient- 
based dynamics which are persistent for some period result in damped oscillations (Figs. 3d 
and 3e). There are cases where the damped oscillations illustrate multipotency (or 
pluripotency depending on the GRN), which is represented by the equal probabilities of 
differentiating towards all the considered cell types (Fig. 3d). In some cases, the interaction 
coefficients vary with different timescales, resulting in partial differentiation and sometimes 
in the reversal of the status of the initially dominant GRF (Fig. 3e). Note that if the 
dampening of the oscillations is fast, the initial oscillations can be unnoticeable yet can still 
activate suppressed genes (Supplementary Fig. 19). 
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Discussion . Various studies have attempted to model the cell differentiation process, but 
there are still more to uncover in epigenetics. Further theoretical prediction and experimental 
validation are needed to fully explain cell-fate determination and reprogramming. Varying 
the efficiency of GRF in expressing a gene (/?), the degradation rate (p), or the constitutive 
growth rate (g) is a straightforward technique in stimulating the activation or deactivation of a 

6 7 10 17 20 

GRF and its corresponding gene ' ' ' ' . However, regulating the repression strength of 
GRFs (y,y) has not been explored, and we have shown numerical illustrations where variations 
in this GRF-GRF repression affect the qualitative behavior of the cell differentiation system 
(Figs. 2 and 3). We are able to replicate Waddington's model using a single set of equations 
with many GRFs involved. A GRN with only two nodes generally cannot describe the 
sequential bifurcation of canals and the oscillations in cell-fate determination. The different 
timescales involved in gene interaction influence the outcome of cellular regulation 5 ' 21 " 24 . 

One of the aims of this study is to spur more discussions on non-equilibrium 
dynamics and oscillations arising from high-dimensional asymmetric systems, which can 
broaden our understanding about the mechanisms of gene regulation. Reversal of the route 
from differentiated state to pluripotency is previously thought to be impossible but now many 
dedifferentiation techniques have been proposed and tested by experiments 9 ' 25 ' 26 . We propose 
another alternative technique for cellular reprogramming, which is by rewiring the GRN to 
have a repressilator-type network, possibly with the aid of external stimulus and stochastic 

1 7 

noise . External stimulus can be introduced to weaken the repression in one loop (Fig. 3a). 
Our numerical predictions can help design cellular engineering strategies for generating 

17 27 28 

induced multipotent stem cells (or pluripotent stem cells depending on the GRN ' ' ) using 
the oscillations that can activate silenced genes. Several studies have discussed that 
oscillating GRF expression is indeed an attribute of progenitor cells 29 " 31 , thus supporting our 
claim. However, note that in reprogramming back to pluripotency, we also need to assure 
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activation of defined factors in the core pluripotency circuitry, such as Oct4, Sox2 and 
Nanog 17 ' 25 . 

The oscillator motifs (e.g., repressilators) which are part of a larger GRN contribute to 
the fluctuations observed in gene regulation dynamics. In fact, there are many types of 
oscillators " . The oscillations generated by each oscillator motif when combined are often 
interpreted as stochastic noise or as chaos. However, note that the combined large and small 
oscillations are not entirely stochastic fluctuations, especially when the detected noise is part 
of the gene regulation system and not just coming from random sources 35 ' 36 . Another 
perspective is that the oscillator motif of a larger GRN can be used for artificial 
transdifferentiation by generating oscillations that can prompt cell-lineage switching, similar 
to what stochastic fluctuations can do 36 ' 37 . In fact, transdifferentiation between related cell 
types branching from one lineage can be more straightforward compared to dedifferentiation 
to pluripotency 38 . 

From our simulations, we formulate some conjectures: (i) The dedifferentiation 
caused by abnormal oscillators (e.g., aberrant repressilator-type network) play a role in the 
existence of cancer stem cells and mutator phenotype 39 " 43 . Abnormal changes in the structure 
of GRF-GRF interaction, such as abnormal timescale factor and abnormal weakening of 
repression links, can lead to disease. Indeed, partially reprogrammed cells and excessive 
plasticity can cause cancer 38 ' 42 ' 44 ' 45 . Moreover, it is also possible that these oscillators play a 
part in epigenomic reprogramming and influence transgenerational epigenetic inheritance 46 . 
Abnormalities in the GRF-GRF interaction could be passed-on to offspring, (ii) We can 
reprogram cells back to pluripotency by regulating the wiring of the GRN. This implies that 
there are no unique reprogramming factors, and we can reprogram cells using any regulatory 

38 

factor as long as it can lower the "gravity" of the epigenetic landscape . In this paper, we 
have shown theoretical predictions of novel developments in the theory of gene regulatory 
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networks. One experimental approach to demonstrate our numerical predictions is by 
employing synthetic biology techniques 19 ' 47 . 

In reality, the temporal transformation in the epigenetic landscape is due to multiple 
intrinsic and extrinsic factors. Here we only consider changes in the GRF-GRF interaction 
coefficient y,y but we should not disregard that gene regulation consists of the interplay among 
many factors and processes. For example, GRF-GRF interaction can be regulated not only 
through yij but also through the modifications in the maximal growth rate ifi+g) or through 

6 7 10 17 ~M ) 

the degradation rate (p) ' ' ' ' . Increasing the maximal growth rate or decreasing the 
degradation rate of a certain GRF enhances the steady state strength of the GRF, which in 
turn intensifies the repression of the other GRFs. Furthermore, the dynamics observed from 
empirical data combine the various effect of many parameters. For example, a decline in the 
strength of a GRF suggests various possible reasons, such as due to an increased degradation 
rate or due to an increased repression by an antagonist GRF. Hence, we need to interpret data 
by considering all possible factors. There are other mechanisms not explicitly discussed in 
this paper, such as spatial pattern of cell differentiation, chromatin remodeling and DNA 
methylation 48 " 50 . 

In summary, our simulations predict the following outcomes: First is the sequential 
branching of lineages in cell-fate determination that portrays differentiation from pluripotent 
state to transient states (lineage progenitors) towards specialized cell types. Second is the 
cellular reprogramming driven by oscillations generated by a repressilator-type network, 
which is a possible technique for dedifferentiation or transdifferentiation. A two-variable 
switch-like model usually cannot illustrate the branching phenomena and sustained 
oscillations, but a high-dimensional model with asymmetric reciprocal interaction between 
GRFs can. Oscillatory behavior cannot be taken for granted because this could explain 
peculiar dynamics related to the epigenetic machinery of organisms, such as dedifferentiation 
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as part of regenerative process. Oscillations are also involved in pattern formation, circadian 

30 33 37 5 1 55 

rhythms, and the progression of diseases ' ' ' " . Network motifs, such as the repressilators, 
that are part of a larger GRN induce functional fluctuations necessary for tissue development 
and cellular engineering. Investigating the dynamics of these network motifs can be helpful in 
drug discovery 56 . 



Methods . In our simulations, we use the following differential equation model (see Box 1): 



dX ; 



fix, 

1 + x i + XrijXj 



dt + cr.dW, i = 1,2, ...,n. 



(1) 



Note that we restrict our simulations to specific parameter values such as P=l, p=0.05 and 
Yij - a, /+ u 2 ■ We suppose all GRFs have the same value of /? and p to highlight the effect of 

time-varying y,y. The term a A dW represents Gaussian white noise with amplitude <j A . Let a A =0 
and oa=0.5 for deterministic and stochastic simulations, respectively. The noise term 
approximates multiple heterogeneous sources of additive random fluctuations. 

Time evolution of the interaction coefficient. The value of y tj - y 1+u 2 is updated using the 
following gradient function: 



du i = exp( —s f t ) 



8U: 



/3X, 



j*i 1 + W; 



2*; 



dt + cr.dW 



(2) 



exp( -Sjt ) 



2A 



tll + { Ui +A) l^i + K- 



Af Xj 



dt + a„dW. 



Equation (2) is used for finding relative optimum and is similar to the trait dynamics 

57 58 

frequently used in evolutionary biology ' . Increasing the value of u, decreases the value of 
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repression coefficient y,y (j^i). Hence, the variable w, can be defined as the time-varying 
attribute for maximizing the strength of the GRF X,. However, note that increasing the value 
of Ui does not always result in an increased equilibrium value of X» especially when the initial 
condition and other parameter values do not allow significant changes in the epigenetic 
landscape in favor of Xj. 

The timescale factor is represented by exp(-Sit) with decline rate e„ as described in 
various studies 5 ' 11 . As time progresses (e.g., as cell matures), the timescale factor declines and 
the value of w, leads to equilibrium. In addition, the dynamic parameter is initialized with 
value Ui(0)=0.001 for all i. Let o u =0 and o u =0.01 for deterministic and stochastic simulations, 
respectively. For simplicity, we approximate the partial derivative in Equation (2) using 
central difference formula with A =0.001. 

Quasi-potential, cell type probability and entropy. The quasi-potential (<P) of the landscape 
is computed as 



d0 
dt 



v dt j 



(3) 



where the dX / dt is deterministic 6 ' 16 . The probability of differentiating to cell type i or the 
expected proportion of cells committed towards cell type i is P. — x yL„ . To visualize the 

canalization in Waddington's illustration, we use three coordinate axes: time, quasi-potential, 

and cell type probability. We also compute for the entropy 59 defined by 

E = -Y d P i log{P i ),P i *0. (4) 

i 

Numerical method. The ordinary and stochastic differential equations are solved using 
Runge-Kutta 4 and Euler-Maruyama with 0.01 as step size, respectively. For supplementary 
mathematical discussions about the differential equation model (Box Eq. 1), refer to related 
literatures 10 ' 17 ' 20 . 
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Notes on abstraction. One of the advantages of the model in Box 1 is that it is 
straightforward, thus any peculiar dynamics can be clearly interpreted. However, the model is 
an abstraction of the cell differentiation process. Hence, we focus on the qualitative dynamics 
rather than on exact quantitative values. Likewise, the qualitative dynamics arising from the 
model are possible to arise in a more complex system that contains the minimal GRN as a 
subnetwork. 
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Figure legends 

Figure 1. A sketch of the epigenetic landscape, and a gene regulatory network (GRN) 
for cell-fate determination, (a) Epigenetic landscape adapted from Waddington's 
illustration 1 . The branching canals depict the various cell lineages towards different fates (cell 
types/phenotypes). The cell fates are illustrated as valleys and traditionally represented as 
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mathematical attractors (see Supplementary Figs. 1 and 2). The differentiating cells, 
illustrated as balls, trail a chosen canal towards a specific valley. The canal is chosen based 
on the landscape's potential (similar to gravitational potential) such that the steeper pathway 
and deeper valley are preferred . The canals are separated by ridges that restrain cells to 
switch lineages. The blue pegs (GRFs) and blue strings (GRF-GRF interaction) alter the 
height of ridges and depth of valleys. The height of ridges and depth of valleys vary through 
time and affect the route of the differentiating cells, (b) A minimal GRN that characterizes 
decision switches in cell-fate determination. We suppose that a node represents a GRF 
involved in expressing a specific cell type, such as master switch genes, transcription factors, 
or coarse-grained modules of a larger GRN that can be simplified as one node. Each node has 
auto-activation as represented by the arrows; while the interaction links among GRFs is of 
repressive behavior (represented by the bars). The repressive behavior of GRFs denotes that a 
mature cell only expresses one phenotype and hinders the expression of multiple phenotypes. 
Strength of repression is not necessarily reciprocal and a one-way repression is possible. The 
auto-activation and repression can be direct or indirect. Examples of GRNs of this type are 

10 22 27 

discussed in various literatures ' ' . 

Figure 2. Illustrations of cells trailing the branching pathways in the epigenetic 
landscape (only the deterministic path is shown; see Supplementary Figs. 3, 5, 6, 8 and 10 
for stochastic simulations and for the initial condition and parameter values used). Note that 
the branch endpoints are the coordinates of an equilibrium point. The different numbers of 
endpoints are due to the different timescale factors used, (a) Five different endpoints. (b) No 
branching and only one endpoint (graphs for GRF 1-5 are superimposed to each other). The 
state with equal cell type probabilities represent multipotent or pluripotent cells 
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(undifferentiated state), (c-e) Two, three and four endpoints, respectively, (f) Time series of 
the entropy levels for the system in Figs. 2a to 2e. The entropy decreases as cells differentiate. 

Figure 3. Oscillating pathways in the epigenetic landscape generated by repressilator- 
type networks, (a) A repressilator-type network with a strong negative feedback loop, n=3 
(this can be extended for any n). The strengths of repression in one loop (solid black bars) are 
stronger than the reverse loop (broken bars). The red bars represent inhibition of repression. 
Note that in our model, only three GRFs are needed to generate sustained oscillations. This 
repressilator-type network has auto-activation unlike the repressilator by Elowitz and 
Leibler 19 . (b-c) Examples of oscillating pathways. There are no deep valleys only continuous 
zig-zag canals (see Supplementary Figs. 13-15 and 17 for the stochastic simulations and for 
the parameter values used), (b) n=5; X5 is initially silenced, (c) n=4; X4 is initially silenced, 
(d) Damped oscillations towards multipotency or pluripotency. The rates of decline of the 
timescale factors are all equal to 0.001. (e) Damped oscillations resulting in partial 
differentiation and reversal of dominant GRF. The rates of decline of the timescale factors are 
not all equal. The initial dominant regulatory factor is GRF 5 but eventually becomes inferior 
as oscillations dampen. 

Supplementary Materials 

Supplementary Figs. 1 to 21 
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Supplementary Figure 1. Example of phase portrait of the ordinary differential 
equation model (Box Eq. 1 in the main text) with more than one stable equilibrium 
point (multistable). The convergence to an equilibrium point depends on the initial 
condition. The coordinates of an equilibrium point (e.g., Ai* and A2*) are the branch 
endpoints of the pathways. 
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Supplementary Figure 2. The changes in the value of interaction coefficient yij drive 
the transformation of the epigenetic landscape affecting the fate of the differentiating 
cell, (a-b) The initial trajectory of the solution to the differential equation model (Box 
Eq. 1 in the main text) converges to a state with equal X1-X2. However, a slight 
modification in the topography of the landscape steers the trajectory towards a state 
with Xi>X2, hence branching arises, (c) Zooming in the phase portrait in Fig. 2b to 
visualize the transformation of the landscape (horizontal axis is Xr, vertical axis is 
X2). Parameters are an- an=l, ei-0.001 and E2-0.005. 
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Supplementary Figure 3. Sample stochastic paths of the system in main text's Fig. 2a 
(five endpoints when deterministic). The initial condition is Xi=l for all i. The 
parameter values are 0jj=l/8 for all i,j, ei=0.0010, £2=0.0060, £3=0.0085, £4=0.0090 and 
e 5 =0.0095. 
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Supplementary Figure 4. Time evolution of u, for the deterministic system in main 
text's Fig. 2a. Timescale factor decline rates are ei-0.0010, £2=0.0060, £3=0.0085, 
£4=0.0090 and £ 5 =0.0095. 
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Supplementary Figure 5. Sample stochastic paths of the system in main text's Fig. 2b 
(one endpoint, undifferentiated state). The initial condition is Xi-1 for all i. The 
parameter values are aij-1/8 and Uij does not evolve for all 
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Supplementary Figure 6. Sample stochastic paths of the system in main text's Fig. 2c 
(two endpoints). The initial condition is X;=2 for all i. The parameter values are 
aij=l/8 for all ei=0.0010, e 2 =0.0012, E3=0.0013, e 4 =0.0100 and e 5 =0.0500. 
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Supplementary Figure 7. Time evolution of u% for the deterministic system in main 
text's Fig. 2c. Timescale factor decline rates are ei-0.0010, £2=0.0012, £3=0.0013, 
£4=0.0100 and £ 5 =0.0500. 




Supplementary Figure 8. Sample stochastic paths of the system in main text's Fig. 2d 
(three endpoints). The initial condition is Xt-1 for all i. The parameter values are 
aij=l/8 for all i,j, ei=0.0010, e 2 =0.0015, e3=0.0090, e 4 =0.0200 and £ 5 =0.0600. 
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Time 

Supplementary Figure 9. Time evolution of u% for the deterministic system in main 
text's Fig. 2d. Timescale factor decline rates are ei-0.0010, £2=0.0015, £3=0.0090, 
£4=0.0200 and £ 5 =0.0600. 
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Supplementary Figure 10. Sample stochastic paths of the system in main text's Fig. 
2e (four endpoints when deterministic). The initial condition is Xi=l for all i. The 
parameter values are for all ei=0.0010, £ 2 =0.0060, £3=0.0090, £4=0.0150 and 

£5=0.0200. 
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Supplementary Figure 11. Time evolution of m for the deterministic system in main 
text's Fig. 2e. Timescale factor decline rates are £i=0.0010, £2=0.0060, £3-0.0090, 
£4=0.0150 and £ 5 =0.0200. 




Supplementary Figure 12. Suppose the parameter values are ai2=a23=a3i=3, ai3=a2i=a32=y, ei-0.001, £2=0.001 and £3=0.001. The initial 
condition is Xi= X2=0 and X3=5. Decreasing the value of y creates a repressilator network such that one repression loop is stronger 
than the reverse loop. This results in an attracting oscillatory behavior that activates suppressed genes. 
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Supplementary Figure 13. Sample stochastic paths of the system in main text's Fig. 
3b, n-5. The initial condition is Xt-0 for i-1,2,3,4 and Xs-15. The parameter values 
are an- a^- 034- #45= asi-5 and aij-0.01 for the other The parameter Wj does not 
evolve for all 
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Supplementary Figure 14. Sample stochastic paths of the system in main text's Fig. 
3c, n-4. The initial condition is Xi-0 for i=l,2,3 and X4-I. The parameter values are 
au= «23= «34= U4i=5, au= «24= «3i= 042=1 and aij=0.01 for the other The parameter un- 
does not evolve for all 
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Supplementary Figure 15. Sample stochastic paths of the system in main text's Fig. 
3d. The initial condition is X;=0 for i=l,2,3,4 and Xs-15. The parameter values are 
an- U23- as4- U45- asi-S, aij-0.01 for the other i,j, ei-0.001, £2=0. 0.001, £3=0. 0.001, 
£4=0.001 and £5=0.001. 
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Supplementary Figure 16. Time evolution of m for the system in main text's Fig. 3d. 
Timescale factor decline rates are Eij-0.001 for all 
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Supplementary Figure 17. Sample stochastic paths of the system in main text's Fig. 
3e. The initial condition is Xt-0 for i=l,2,3,4 and X5-I5. The parameter values are au- 
ai3= 034= a 4 5= 0151=5, cnj=0.01 for the other ei=0.001, £ 2 =0. 0.001, £3=0. 0.001, £4=0.010 
and £5=0.005. 
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Supplementary Figure 18. Time evolution of w for the system in main text's Fig. 3e. 
Timescale factor decline rates are ei-0.001, £2=0. 0.001, £3-0. 0.001, £4-0.010 and 
£5=0.005. 
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Supplementary Figure 19. Sample stochastic paths of the system with slow 
timescale factor decline rate, a-0.0001 for all i. Here dampening of the oscillations is 
fast and the initial oscillations are unnoticeable; however, repressilator-type network 
is still able to activate suppressed genes. The initial condition is X;=0 for i=l,2,3,4 and 
X5=15. The parameter values are an- «23= 034= 045= asi-l and aij-0.1 for the other 
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Supplementary Figure 20. Time series of entropy levels for the systems in the main 
text's Fig. 3d and 3e. 
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Supplementary Figure 21. Deterministic system with 10 GRFs. The initial condition 
is Xi-1 for all i. The parameter values are cnj=l/8 for all ei=0.0001, £2-0.0021, 
£3=0.0042, £4=0.0056, £ 5 =0.0068, £(,=0.0077, £7=0.0081, £ 8 =0.0083, £ 9 =0.0085 and 
£w=0.0089. 



